
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

  ncfnama = "radon_mine_E0_1999-2003_dm_ymm.nc" 
  ncfnamb = "radon_mine_E1_1999-2003_dm_ymm.nc" 
  ncfnamc = "radon_mine_E2_1999-2003_dm_ymm.nc" 
  ncfnamd = "radon_mine_E3_1999-2003_dm_ymm.nc" 
  ncfname = "radon_obs_mm.nc" 
  ascnama = "radon_siteplot.dat" 
  ascnamb = "radon_sitename.dat" 
  ascnamc = "radon_sitebad.dat" 

  fla = addfile(ncfnama,"r") 
  flb = addfile(ncfnamb,"r") 
  flc = addfile(ncfnamc,"r") 
  fld = addfile(ncfnamd,"r") 
  fle = addfile(ncfname,"r") 

  rma = fla->radon
  rmb = flb->radon
  rmc = flc->radon
  rmd = fld->radon
  rme = fle->radon

  kplot = asciiread(ascnama,-1,"integer") 
  sname = asciiread(ascnamb,-1,"string") 
  kbad  = asciiread(ascnamc,-1,"integer") 

  ns = dimsizes(kplot) 

  rma!1 = "site" 
  rmb!1 = "site" 
  rmc!1 = "site" 
  rmd!1 = "site" 

  ;do is = 0,ns-1 
  ;   if(kbad(is).eq.0) then 
  ;      rma(:,is) = -999  
  ;      rmb(:,is) = -999  
  ;      rmc(:,is) = -999  
  ;      rmd(:,is) = -999  
  ;      rme(:,is) = -999  
  ;   end if  
  ;end do 
 
  sma = dim_avg(rma(site|:,time|:))  
  smb = dim_avg(rmb(site|:,time|:))  
  smc = dim_avg(rmc(site|:,time|:))  
  smd = dim_avg(rmd(site|:,time|:))  
  sme = dim_avg(rme(site|:,time|:))  

  nt = 12 

  ;do is = 0,ns-1 
  ;;if(kplot(is).eq.0) then 
  ;print("-") 
  ;print("obs: "+sname(is)+"  "+rme(:,is)+"   "+rma(:,is)+"   "+rmd(:,is)) 
  ;print("-") 
  ;if(.not.ismissing(sme(is)) .and. .not. ismissing(-rme(0,is)) ) then 
  ;;do it = 0,nt-1 
  ;if(abs(sme(is)-rme(0,is)).le.1.) then 
  ;   print("is: "+sname(is)) 
  ;   rme(:,is) = -999 
  ;   rma(:,is) = -999 
  ;   rmb(:,is) = -999 
  ;   rmc(:,is) = -999 
  ;   rmd(:,is) = -999 
  ;end if 
  ;;end do 
  ;end if 
  ;;end if 
  ;end do 

  sma = dim_avg(rma(site|:,time|:))  
  smb = dim_avg(rmb(site|:,time|:))  
  smc = dim_avg(rmc(site|:,time|:))  
  smd = dim_avg(rmd(site|:,time|:))  
  sme = dim_avg(rme(site|:,time|:))  

  ;print("radon: "+sma+"   "+sname) 

  is = 0 

  do is = 0,ns-1  

  print("kplot,is : "+kplot(is)+"  "+is) 

  ;if( kbad(is).ne.0) then
  
  print(""+sname(is)) 
 
  L_log = 0
 
  plot = new( 2,graphic )

  LeftString  = (/"Radon concentration (mBq m~S~-3~N~) "+sname(is)/) 

  XAxisString = "Observation"
  YAxisString = "Model"
    
  outname = "mm_scatter"  

  wks = gsn_open_wks( "ps",outname+sname(is) )

   x1 = ndtooned(rme(:,is)) 
   y1 = ndtooned(rma(:,is)) 
   y2 = ndtooned(rmd(:,is)) 
   y3 = ndtooned(rmc(:,is)) 
   ;x1 = ndtooned(rme(:,:)) 
   ;y1 = ndtooned(rma(:,:)) 
   ;y2 = ndtooned(rmd(:,:)) 
   ;y3 = ndtooned(rmc(:,:)) 

   x1 = where(x1.lt.0.,-9999999.,x1)
   y1 = where(y1.lt.0.,-9999999.,y1)
   y2 = where(y2.lt.0.,-9999999.,y2)
   y3 = where(y3.lt.0.,-9999999.,y3)

   ;print("x1,y1,y2: "+x1+"   "+y1+"   "+y2) 

   ic = 0
   ia = 0
   
   do i = 0, dimsizes(x1)-1
   if(.not.ismissing(x1(i)).and..not.ismissing(y1(i))) then  
   if(x1(i) .ge. 0.00001) then 
      if(y1(i)/x1(i).le. 0.5 .or. y1(i)/x1(i).ge. 2.0 ) then 
      	 ic = ic + 1
      end if
      ia = ia + 1
   end if
   end if
   end do

   if(ia.ne.0) then 
   fy1 = 100. - int2flt(ic) / int2flt(ia) * 100.
   else
   fy1 = 0. 
   end if 

   
   ic = 0
   ia = 0
   
   do i = 0, dimsizes(x1)-1
   if(.not.ismissing(x1(i)).and..not.ismissing(y1(i))) then  
   if(x1(i) .ge. 0.00001) then 
      if(y2(i)/x1(i).le. 0.5 .or. y2(i)/x1(i).ge. 2.0 ) then 
      	 ic = ic + 1
      end if
      ia = ia + 1
   end if
   end if
   end do 
   
   if(ia.ne.0) then 
   fy2 = 100. - int2flt(ic) / int2flt(ia) * 100.
   else
   fy2 = 0. 
   end if 
 
;; fy3 
   
   ic = 0
   ia = 0

   do i = 0, dimsizes(x1)-1
   if(.not.ismissing(x1(i)).and..not.ismissing(y1(i))) then
   if(x1(i) .ge. 0.00001) then
      if(y3(i)/x1(i).le. 0.5 .or. y3(i)/x1(i).ge. 2.0 ) then
         ic = ic + 1
      end if
      ia = ia + 1
   end if
   end if
   end do 

   if(ia.ne.0) then
   fy3 = 100. - int2flt(ic) / int2flt(ia) * 100.
   else
   fy3 = 0. 
   end if 


   print(""+fy1+" "+fy2+" "+fy3)
  

   x1@_FillValue = -9999999.
   y1@_FillValue = -9999999.
   y2@_FillValue = -9999999.
   y3@_FillValue = -9999999.
;;print(x1)

   fr1 = esccr(x1,y1,0)
   fr2 = esccr(x1,y2,0)
   fr3 = esccr(x1,y3,0)
 

   x1 = where(x1.lt.0.,9999999.,x1)
   y1 = where(y1.lt.0.,9999999.,y1)
   y2 = where(y2.lt.0.,9999999.,y2)
   y3 = where(y3.lt.0.,9999999.,y3)


   
   ResL          = True
   ResL@gsnDraw  = False
   ResL@gsnFrame = False

   ResL@tiYAxisFontHeightF     = 0.016
   ResL@tmYLLabelFontHeightF   = 0.016
   ResL@tmXBLabelFontHeightF   = 0.016
   ResL@gsnLeftStringFontHeightF   = 0.016
   ResL@gsnRightStringFontHeightF  = 0.016
 
 
   ResL@tiXAxisFontHeightF     = 0.016
   ResL@tiYAxisFontHeightF     = 0.016
   ResL@tiXAxisString       = XAxisString
   ResL@tiYAxisString       = YAxisString
   
   if(L_log.eq.0) 
   ResL@trXLog              = True
   ResL@trYLog              = True
   ResL@trXMinF             = 50
   ResL@trXMaxF             = 50000
   ResL@trYMinF             = 50
   ResL@trYMaxF             = 50000
   else
   ResL@trXLog              = False
   ResL@trYLog              = False
   ResL@trXMinF             = 0.
   ResL@trXMaxF             = 5000
   ResL@trYMinF             = 0.
   ResL@trYMaxF             = 5000
   end if
   	

   ResL@gsnLeftString  = LeftString
   
   ResL@trXReverse          = False

   ResL@xyMarkLineModes   = "Markers"             ; choose which have markers
   ResL@xyMarkers         = (/1/)                 ; choose type of marker  
   ResL@xyMarkerSizeF     = 0.025                  ; Marker size (default 0.01)
   ResL@xyMarkerColor     = (/"red"/)             ; Marker color



   xx = new((/dimsizes(x1),2/),"float") 
   yy = new((/3,dimsizes(y1)/),"float") 
  
   ;;xx(:,0) = x1(:)  
   ;;xx(:,1) = x1(:)  

   yy(0,:) = y1(:) 
   yy(1,:) = y2(:) 
   yy(2,:) = y3(:) 
 
   plot1 = gsn_csm_xy(wks,x1,y1,ResL)



   txres               = True                     ; text mods desired
   txres@txFontHeightF = 0.015                     ; font smaller. default big
   txres@txJust        = "CenterLeft"              ; Default is "CenterCenter".

   gsn_text_ndc(wks,"~F29~P2: " + sprintf("%3.2f",fy1) + " CR: " + sprintf("%3.2f",fr1),0.6,0.20,txres) 
   gsn_text_ndc(wks,"~F29~P2: " + sprintf("%3.2f",fy2) + " CR: " + sprintf("%3.2f",fr2),0.6,0.23,txres) 
   gsn_text_ndc(wks,"~F29~P2: " + sprintf("%3.2f",fy3) + " CR: " + sprintf("%3.2f",fr3),0.6,0.27,txres) 
;; gsn_text_ndc(wks,"~F29~CR: " + sprintf("%3.2f",fr1),0.7,0.20,txres) 
;; gsn_text_ndc(wks,"~F29~P2: " + sprintf("%3.2f",fy1),0.7,0.20,txres) 
;; gsn_text_ndc(wks,"~F29~P2: " + sprintf("%3.2f",fy2),0.7,0.30,txres) 
;; gsn_text_ndc(wks,"~F29~CR: " + sprintf("%3.2f",fr2),0.7,0.35,txres) 

   ;;drawNDCGrid(wks) 
 
   ResL@xyMarkerColor     = (/"blue"/)          ; Marker color
   plot2 = gsn_csm_xy(wks,x1,y2,ResL)

ResL@xyMarkerColor     = (/"green"/)          ; Marker color
   plot3 = gsn_csm_xy(wks,x1,y3,ResL)
 
   overlay(plot1,plot2)
   overlay(plot1,plot3)
   
   resl = True

   resl@gsLineColor = "black"
   resl@gsLineLabelString= ""

   xpts = (/1.,1./)
   ypts = (/1.,1./)


   ypts(0) = ResL@trYMinF
   ypts(1) = ResL@trYMaxF
   xpts(0) = ResL@trXMinF
   xpts(1) = ResL@trXMaxF

   resl@gsLineThicknessF = 0.1
   resl@gsLineDashPattern = 0

   duna=gsn_add_polyline(wks,plot1,xpts(0:1),ypts(0:1),resl)


   resl@gsLineThicknessF = 1.0
   resl@gsLineDashPattern = 11

   ypts(0) = ResL@trYMinF
   ypts(1) = 0.5 * ResL@trYMaxF
   xpts(0) = 2. * ResL@trXMinF
   xpts(1) = ResL@trXMaxF

   dunb=gsn_add_polyline(wks,plot1,xpts(0:1),ypts(0:1),resl)

   ypts(0) = 2. * ResL@trYMinF
   ypts(1) = ResL@trYMaxF
   xpts(0) = ResL@trXMinF
   xpts(1) = 0.5 * ResL@trXMaxF

   dunc=gsn_add_polyline(wks,plot1,xpts(0:1),ypts(0:1),resl)


 
   ResP                            = True
;; ResP@gsnMaximize                = True
;; ResP@gsnPaperMargin             = 0.2
;; ResP@gsnPanelXWhiteSpacePercent = 5
;; ResP@gsnPanelYWhiteSpacePercent = 5

   gsn_panel (wks,plot1, (/1,1/),ResP)
     
   frame(wks)

   ;else
;
;   end if 

   end do 



